#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 24 22:08:30 2021

@author: liqingsimac
"""

'''
########################################################
##Example-8-1-1
import numpy as np
import matplotlib.pyplot as plt

def myfun(x,y):
    numerator = x+y*(x**2+y**2-1)
    denominator = -y+x*(x**2+y**2-1)
    if denominator !=0:
        return numerator/denominator
    else:
        return np.inf

b=1.5; N=17; 
stk=0.1 # length of each stick (line element)

x=np.linspace(-b,b,num=N)
y=np.linspace(-b,b,num=N)

fig=plt.figure()
ax=fig.add_subplot(111)
for k in range(N):
    for m in range(N):
        xp=x[k]
        yp=y[m]
        slope=myfun(xp,yp)
        if (slope==np.inf):
            ax.plot([xp,xp],[yp-stk/2,yp+stk/2],color='blue')
        else:
            dx=stk/np.sqrt(1+slope**2)
            dy=dx*slope
            ax.plot([xp-dx/2,xp+dx/2],[yp-dy/2,yp+dy/2],color='blue')

t=np.linspace(0,2*np.pi,101)
x1=np.cos(t)
y1=np.sin(t)
ax.plot(x1,y1,'c--')

ax.plot(0,0,'g.',markersize=10)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
#fig.savefig('ode-example-8-1-1.png')
########################################################
'''

'''
########################################################
##Example 8.3.1
import numpy as np
import matplotlib.pyplot as plt

def myfun(C1,C2,t):
    x=3*C1*np.exp(3*t)+C2*np.exp(-4*t)
    y=C1*np.exp(3*t)-2*C2*np.exp(-4*t)
    return x,y

m=20
a=0.3; N=101
t=np.linspace(-a,a,N)
C1=[-4,-2,0,2,4]
C2=C1

fig=plt.figure()
ax=fig.add_subplot(111)

for k in range(len(C1)):
    for m in range(len(C2)):
        x,y=myfun(C1[k],C2[m],t)
        ax.plot(x,y,'b-')
ax.set_aspect('equal')
ax.set_xticks([])
ax.set_yticks([])
#ax.set_xlim(-m,m)
#ax.set_ylim(-m,m)
#fig.savefig('ode-example-8-3-1.png')
########################################################
'''

'''
########################################################
##Example 8.3.2
import numpy as np
import matplotlib.pyplot as plt

def myfun(C1,C2,t):
    x=C1*np.exp(3*t)
    y=C1*np.exp(3*t)/2+C2*np.exp(t)
    return x,y

a=0.3; N=101
t=np.linspace(-8*a,a,N)
C1=[-8,-4,0,4,8]
C2=[-8,-4,0,4,8]

fig=plt.figure()
ax=fig.add_subplot(111)

for k in range(len(C1)):
    for m in range(len(C2)):
        x,y=myfun(C1[k],C2[m],t)
        ax.plot(x,y,'b-')
ax.set_aspect('equal')
ax.set_xlim([-10,10])
ax.set_ylim([-10,10])
ax.set_xticks([])
ax.set_yticks([])
#ax.set_yticks([-10,-5,0,5,10])
#fig.savefig('ode-example-8-3-2.png')
########################################################
########################################################
########################################################
########################################################
'''

'''
########################################################
##8.3：图8-6-7：星形结点
import numpy as np
import matplotlib.pyplot as plt

m=5
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

#A=[-1,0,0,-1] #星形结点、稳定
A=[1,0,0,1] #星形结点、不稳定
[a,b,c,d]=A

U=a*X+b*Y
V=c*X+d*Y
C=(X**2+Y**2)/np.max(X**2+Y**2)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
t=np.linspace(-3,3,21)
C1=[-4,-2,-1,0,2,5]
C2=[-4,-2,-1,0,2,5]
for k1 in range(len(C1)):
    for k2 in range(len(C2)):
        ax.plot(C1[k1]*np.exp(a*t),C2[k2]*np.exp(d*t),'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
'''

'''
########################################################
##8.3：图8-8-9：两向结点
import numpy as np
import matplotlib.pyplot as plt

m=5
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

#A=[-1,0,0,-2] #两向结点、稳定
A=[2,0,0,1] #两向结点、不稳定
[a,b,c,d]=A

U=a*X+b*Y
V=c*X+d*Y
C=(X**2+Y**2)/np.max(X**2+Y**2)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
t=np.linspace(-3,3,21)
C1=[-4,-2,-1,0,2,5]
C2=[-4,-2,-1,0,2,5]
for k1 in range(len(C1)):
    for k2 in range(len(C2)):
        ax.plot(C1[k1]*np.exp(a*t),C2[k2]*np.exp(d*t),'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
#'''

'''
########################################################
##8.3：图8-6：鞍点
import numpy as np
import matplotlib.pyplot as plt

m=5
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

#A=[-1,0,0,1] # 鞍点
#A=[1,0,0,-1] # 鞍点
[a,b,c,d]=A

U=a*X+b*Y
V=c*X+d*Y
C=(X**2+Y**2)/np.max(X**2+Y**2)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
t=np.linspace(-3,1,21)
C1=[-4,-2,-1,0,2,5]
C2=[-4,-2,-1,0,2,5]
for k1 in range(len(C1)):
    for k2 in range(len(C2)):
        ax.plot(C1[k1]*np.exp(a*t),C2[k2]*np.exp(d*t),'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
'''

'''
########################################################
##8.3：图8-11-12：单向结点
import numpy as np
import matplotlib.pyplot as plt


m=5
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

A=[-1,0,1,-1]
#A=[1,0,1,1]
[a,b,c,d]=A

U=a*X+b*Y
V=c*X+d*Y
C=(X**2+Y**2)/np.max(X**2+Y**2)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
C=[-4,-3,-2,-1,0,1,2,3,4,5]
xg=np.linspace(-m,m,30)
for k in range(len(C)):
    yg=C[k]*xg+xg*np.log(np.abs(xg))/a
    ax.plot(xg,yg,'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
'''

'''
########################################################
##8.3：图8-13-14-15：焦点-中心
import numpy as np
import matplotlib.pyplot as plt

m=10
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

A=[1,-2,2,1] #不稳定焦点
#A=[-1,-2,2,-1] #稳定焦点
#A=[0,-2,2,0] #中心
[a,b,c,d]=A

U=a*X+b*Y
V=c*X+d*Y
C=(X**2+Y**2)/np.max(X**2+Y**2)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
C=[0.1,0.3,0.5,1,2,5,10,15,30] #焦点
#C=[0.5,1,2,3,5,7] #中心
theta = np.linspace(-np.pi*4,np.pi*4,101)
for k in range(len(C)):
    r=C[k]*np.exp(a/c*theta)
    xg=r*np.cos(theta)
    yg=r*np.sin(theta)
    ax.plot(xg,yg,'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
'''

'''
########################################################
##8.1.1：李雅普诺夫稳定性
import numpy as np
import matplotlib.pyplot as plt

m=1.5
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

U=-Y+X*(X**2+Y**2-1)
V=X+Y*(X**2+Y**2-1)
C=(X**2+Y**2)/np.max(X**2+Y**2)

U=(np.abs(U)**(1/3))*np.sign(U)
V=(np.abs(V)**(1/3))*np.sign(V)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
#C=[0.5,1,2,3,5,7]
theta = np.linspace(0,np.pi*2,31)
r=1
xg=r*np.cos(theta)
yg=r*np.sin(theta)
ax.plot(xg,yg,'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
'''

'''
########################################################
## Examle-8-3-1-b：鞍点
import numpy as np
import matplotlib.pyplot as plt

m=15
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

A=[2,3,2,-3] # 鞍点
[a,b,c,d]=A

U=a*X+b*Y
V=c*X+d*Y
C=(X**2+Y**2)/np.max(X**2+Y**2)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
t=np.linspace(-3,3,101)
C1=[-4,-2,-1,0,2,5]
C2=[-4,-2,-1,0,2,5]
for k1 in range(len(C1)):
    for k2 in range(len(C2)):
        xg=3*C1[k1]*np.exp(3*t)+C2[k2]*np.exp(-4*t)
        yg=C1[k1]*np.exp(3*t)-2*C2[k2]*np.exp(-4*t)
        ax.plot(xg,yg,'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
'''

'''
########################################################
## Examle-8-3-2-b：
import numpy as np
import matplotlib.pyplot as plt

m=15
x=np.linspace(-m,m,15)
y=np.linspace(-m,m,15)
[X,Y] = np.meshgrid(x,y)

A=[3,0,2,1] # 单向结点
[a,b,c,d]=A

U=a*X+b*Y
V=c*X+d*Y
C=(X**2+Y**2)/np.max(X**2+Y**2)

fig, ax = plt.subplots()
ax.quiver(X,Y,U,V,C)
ax.set_xticks([])
ax.set_yticks([])
ax.set_aspect('equal')
       
t=np.linspace(-3,3,101)
C1=[-4,-2,-1,0,2,5]
C2=[-4,-2,-1,0,2,5]
for k1 in range(len(C1)):
    for k2 in range(len(C2)):
        xg=C1[k1]*np.exp(3*t)
        yg=C1[k1]*np.exp(3*t)/2+C2[k2]*np.exp(t)
        ax.plot(xg,yg,'-')
ax.set_xlim(-m,m)
ax.set_ylim(-m,m)
########################################################
'''


